Transmission Spectroscopy of the Lowest-density Gas Giant: Metals and a Potential Extended Outflow in HAT-P-67b

In order to isolate the atmospheric signal of the exoplanet, we first needed to remove the telluric and stellar features from the data.

3.1.1. Telluric Correction

3.1.2. Correction of the Rossiter-McLaughlin and Center-to-limb Variation Effects

Inhomogeneities in the spectrum of a star across its disk can significantly impact transmission spectroscopy studies. Of particular relevance to studies at high resolution are the Rossiter-McLaughlin (RM) and center-to-limb variation (CLV) effects. The RM effect occurs when a planet blocks regions of the rotating stellar disk that have different line-of-sight velocities. The CLV effect is the result of changes in the stellar line profiles and intensities between the center and the limb of the stellar disk (Yan et al. 2017). The contamination from the RM and CLV effects is exacerbated in systems like HAT-P-67, where the radial velocity of the planet during transit is close to the line-of-sight velocity of the region of the star occulted by the planet (see, e.g., Casasayas-Barris et al. 2022; Dethier & Bourrier 2023).

We modeled and removed the contribution from the RM and CLV effects to each of our spectra. As in similar works (e.g., Yan et al. 2017; Casasayas-Barris et al. 2019), we produced a grid of synthetic stellar spectra at different limb angles using Spectroscopy Made Easy version 522 (Piskunov & Valenti 2017), the Kurucz ATLAS9 model atmospheres (Castelli & Kurucz 2003), and the Vienna Atomic Line Database (Ryabchikova et al. 2015). We discretized the stellar disk into a series of squared cells (in our case, of side length 0.01 R_⋆), and we assigned to each cell a spectrum interpolated from our grid based on the limb angle of such a cell. We also Doppler shifted each spectrum according to the line-of-sight velocity of the cell. At the orbital phase of each exposure, we integrated the spectra from the cells blocked by the planetary disk, and divided the resulting spectrum by the full-disk spectrum. We continuum normalized the resulting spectra using a median filter, and we broadened them according to the resolution of the corresponding CARMENES channel. This procedure generates a model of the RM and CLV contributions for each exposure, which we then divided out from the data. In Figure 1, we illustrate the effect that the removal of the RM and CLV effects has on the cross-correlation functions (CCFs) of Na i (see Section 3.3 for a description of how we calculate the CCFs).

Zoom In Zoom Out Reset image size

3.1.3. Removal of the Stellar Spectra

3.1.4.  SYSREM as an Alternative Detrending Approach

Following the methodology presented in Bello-Arufe et al. (2022b), we constructed a model of the transmission spectrum for each of the species we were interested in. We assumed an atmosphere of solar metallicity and in chemical equilibrium, and an isothermal profile with a temperature equal to the equilibrium temperature of the planet, T_eq = 1903 K. For each gas, we calculated its concentration with FastChem (Stock et al. 2018), and we computed its absorption cross section with HELIOS-K (Grimm et al. 2021) using the line lists from Kurucz (2018) and assuming a Voigt profile for the absorption lines. Our models also include H⁻ bound-free and free–free absorption (John 1988). At each wavelength, we calculated the transit depth using the formalism in Gaidos et al. (2017) and Bower et al. (2019).

We processed the models similarly to the data in order to turn them into templates for the cross-correlation and likelihood analyses. First, we broadened each model according to the instrumental resolution of each CARMENES channel. Then, we removed from each model its continuum, found using a Gaussian filter with a standard deviation equivalent to that of the Gaussian filter applied to the data.

Our survey focused on atomic and molecular species with a substantial number of detectable lines in the transmission spectrum of a planet like HAT-P-67b and in the CARMENES wavelength range. Our list included Al i, Ca i, Ca ii, CO, Cr i, Fe i, H₂O, K i, Li i, Mg i, Mn i, Na i, Ni i, O i, OH, Sc i, Sc ii, Si i, Ti i, TiO, V i, and VO (Rothman et al. 2010; McKemmish et al. 2016; Kurucz 2018; Polyansky et al. 2018; McKemmish et al. 2019).

Cross-correlation of high-resolution spectra against model templates has become one of the most successful techniques to robustly identify species in exoplanetary atmospheres. This technique has demonstrated its success in the detection of atomic and molecular species, in both transmission and emission spectroscopy (e.g., Snellen et al. 2010; Brogi et al. 2012; Hoeijmakers et al. 2019; Pino et al. 2020; Kesseli et al. 2022).

For each species, we calculated the corresponding CCFs using the expression in Gibson et al. (2020):

$$\begin{eqnarray}&&{\rm{CCF}}=\sum \displaystyle \frac{{f}_{i}{m}_{i}}{{\sigma }_{i}^{2}},\end{eqnarray} \tag{ 1 }$$

where f_i and m_i refer to the values of the spectrum and model template at pixel i, and ${\sigma }_{i}^{2}$ is the variance with time of the values in that pixel. The sum spans all pixels in the corresponding CARMENES channel. We Doppler shifted each model template by velocity offsets Δv_sys ranging from −500 to 500 km s⁻¹ in steps of 1 km s⁻¹. We also scaled the model template by a value between 0 and 1 according to the transit light curve, calculated using batman (Kreidberg 2015).

We then generated the cross-correlation K_p –Δv_sys maps. To do so, we Doppler shifted each CCF by a velocity v_p :

$$\begin{eqnarray}&&{v}_{p}\left(t\right)={K}_{p}\sin \left(2\pi \phi (t)\right),\end{eqnarray} \tag{ 2 }$$

where ϕ(t) is the orbital phase at time t, and K_p is the planet's radial velocity semiamplitude, which we allowed to take values between 0 and 500 km s⁻¹ in steps of 1 km s⁻¹. For each K_p , we summed the Doppler-shifted CCFs and stacked them to generate the K_p –Δv_sys maps. We converted the values in these maps to S/Ns by normalizing them by the standard deviation of those values located at ∣Δv_sys∣ ≥ 50 km s⁻¹, far from the exoplanet signal. Any detection should manifest as a peak near K_p ∼ 147 km s⁻¹ (the planet's radial velocity semiamplitude according to Zhou et al. 2017) and v_sys = 0 km s⁻¹ (as we shifted all spectra to the rest frame of the HAT-P-67 system).

While classical cross-correlation techniques are highly effective to detect atomic and molecular species, their ability for model comparison is more limited. To overcome this issue, recent works have introduced novel Bayesian retrieval frameworks to extract physical parameters from high-resolution ground-based observations (e.g., Brogi & Line 2019; Gibson et al. 2020).

In this work, we adopted the statistical framework presented by Gibson et al. (2020) to constrain the physical extent of the model atmosphere, in addition to K_p and Δv_sys. We computed the log-likelihood as

$$\begin{eqnarray}&&\mathrm{ln}L=\displaystyle \frac{-N}{2}\mathrm{ln}\displaystyle \frac{{\chi }^{2}}{N},\end{eqnarray} \tag{ 3 }$$

where N is the number of data points, and χ² is given by:

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i=1}^{N}\displaystyle \frac{{f}_{i}^{2}}{{\sigma }_{i}^{2}}+{\alpha }^{2}\displaystyle \sum _{i=1}^{N}\displaystyle \frac{{m}_{i}^{2}}{{\sigma }_{i}^{2}}-2\alpha \displaystyle \sum _{i=1}^{N}\displaystyle \frac{{f}_{i}{m}_{i}}{{\sigma }_{i}^{2}}.\end{eqnarray} \tag{ 4 }$$

Here, α is a parameter that provides information on how well our model captures the physical extent of the atmosphere. This parameter is particularly useful to interpret data from highly irradiated gas giants, as these planets often possess extended atmospheres that models struggle to predict (e.g., Hoeijmakers et al. 2019). The final term of this equation is identical to the CCFs from Equation (1). Therefore, as presented in Nugroho et al. (2020), we calculated the log-likelihood as a function of three parameters (α, K_p , and Δv_sys), which led to a three-dimensional data cube for each species, with these three parameters as the axes. We gave α values from 0.1–10, and we kept the same range of values for K_p and Δv_sys as in the previous section. We subtracted the maximum log-likelihood from each data cube and took the exponential to produce a scaled likelihood, with values ranging from 0–1.

Here, we describe the steps we followed to search for hydrogen and helium escape in Hα and the metastable He triplet. These steps were similar to the ones we followed to prepare the spectra for cross correlation, except for some notable differences in the detrending process.

We removed telluric absorption lines using molecfit, as described in Section 3.1.1. The spectral regions around Hα and the He triplet are mostly affected by water telluric absorption. The removal of these absorption lines for a sample spectrum is shown in the first row of Figure 2. We found that molecfit achieved a high-quality modeling of these telluric lines, and hence we did not need to apply the more aggressive SYSREM algorithm.

Zoom In Zoom Out Reset image size

While the region near Hα is free from any significant telluric emission, the sky-monitoring fiber (i.e., fiber B) revealed strong emission lines around the He metastable triplet associated with OH airglow (Oliva et al. 2015). In particular, there are two Λ-split OH doublets from the Q branch that surround the He triplet. The first OH doublet has vacuum wavelengths of 10832.103 and 10832.412 Å. The second OH doublet, about an order of magnitude stronger than the first one, has wavelengths of 10834.241 and 10834.338 Å. The two individual components of this second doublet are not resolved by CARMENES and therefore appear as one strong line in the spectra. The top-right panel of Figure 2 shows the removal of these three telluric emission features from a sample spectrum.

We corrected the OH telluric emission lines using data from fiber B and following a similar modeling approach to Palle et al. (2020a). First, we shifted the spectra of the sky-monitoring fiber to a common wavelength grid through a quadratic interpolation, and we combined them into a median sky spectrum. We then fit three Gaussian functions to the three resolved OH lines, keeping the amplitudes and widths of the two components of the resolved doublet equal to each other. We scaled this model to match the amplitude of the strongest OH feature in each individual sky spectrum, allowing for a small shift in the position of the peaks. Because of the difference in throughput between both fibers, we multiplied the individual sky models by a scaling factor of 0.881 prior to subtracting them from the science spectra. We found that this approach to remove the OH telluric lines provided slightly cleaner spectra in this region than just directly subtracting the sky spectra as described in Section 3.1.1.

We then removed a model of the RM and CLV contributions following Section 3.1.2 and normalized the resulting spectra. As in Palle et al. (2020a), we normalized the data by dividing each individual spectrum in the order that contained the He triplet by the median value in a region of the same order but near the continuum, far from any stellar lines. We repeated this step in the order that contained Hα. We shifted the normalized spectra from the rest frame of the observatory to that of the HAT-P-67 system, as we did in the cross-correlation analysis (Section 3.1.3). After shifting all spectra, we interpolated them onto a common wavelength grid. The second row of Figure 2 shows the normalized spectra in the rest frame of the HAT-P-67 system.

We computed the transmission spectra by dividing the spectra by the master stellar spectrum. As we might be probing the escaping atmosphere, any potential H and He signal may not necessarily follow the planetary orbital velocity and might be broader in wavelength than those originating from the bounded atmosphere. For that reason, we did not include any in-transit spectra in this calculation of the master stellar spectrum. Additionally, as will be evidenced by the transmission light curves later in the analysis, there may be planetary absorption prior to the point of first contact. Therefore, we calculated the master stellar spectrum by averaging only the spectra taken after transit. We show the resulting transmission spectra in the third row of Figure 2.

Finally, we obtained the master transmission spectrum of HAT-P-67b by taking the weighted average of all of the transmission spectra during transit in the planetary rest frame. The last row of Figure 2 shows the transmission spectra stacked in the planetary rest frame.

The cross-correlation analysis reveals the presence of Ca ii and Na i in the atmosphere of HAT-P-67b. Figure 3 presents the cross-correlation K_p –Δv_sys maps for each of these species. The S/N of these detections, as measured at the peak of the signals, is 10.9 in the case of Ca ii, and 4.1 in the case of Na i (see Table 1).

Zoom In Zoom Out Reset image size

Switching from cross-correlation to likelihood maps, which are also shown in Figure 3, allows us to place constraints on K_p , Δv_sys, and α for each species. We present these results in Table 1. As in Nugroho et al. (2020), we calculated these parameters by fitting a Gaussian function to the slice that goes through the maximum scaled likelihood of the data cubes described in Section 3.4. The standard deviation of the Gaussian function gives us an estimate of the uncertainty of the parameter.

The retrieved K_p and Δv_sys values of the Ca ii and Na i detections are consistent with the expected orbital values. The planetary radial velocity semiamplitude (K_p ≈ 147 km s⁻¹, Zhou et al. 2017) is within 1σ of the retrieved K_p values for both detections. Unlike in many other hot and ultrahot gas giants (e.g., Snellen et al. 2010; Casasayas-Barris et al. 2019; Nugroho et al. 2020; Bello-Arufe et al. 2022b; Kesseli et al. 2022; Prinoth et al. 2022; Zhang et al. 2022), we do not measure the blueshifted absorption often associated with day-to-night winds.

The values we retrieve for α inform us of the extent of the atmosphere associated with each species. As shown in Table 1, we find that the retrieved Ca ii signal has the largest value of α; our model underestimates the line contrast of this species by a factor of ∼6. This is similar to results in high-resolution studies of other highly irradiated planets (e.g., Yan et al. 2019; Nugroho et al. 2020; Borsa et al. 2021; Tabernero et al. 2021; Bello-Arufe et al. 2022b), where ionized calcium is often found at very high altitudes, potentially due to photoionization or the presence of an escaping envelope.

Meanwhile, our models overestimate the line contrast of Na i by a factor of 2. In our analysis, we broaden the model templates according to the instrumental resolution, but we ignore other line-broadening mechanisms, such as orbital smearing and rotational broadening, which can bias the retrieved value of α. In each exposure, the planet changes radial velocity by only ∼1 km s⁻¹, so we do not expect orbital smearing to significantly alter α. However, planetary rotation can have a measurable effect on the line depths. For example, rotational broadening of WASP-76b can reduce its line depths by a factor of ∼2 (Casasayas-Barris et al. 2021). While rotational broadening in HAT-P-67b (v_rot ∼ 2 km s⁻¹, assuming tidal locking) is less severe than in WASP-76b (v_rot ∼ 5 km s⁻¹), this mechanism may partially explain why the Na i lines of HAT-P-67b are shallower than we expected. Other factors can further reduce line contrasts. Silicate clouds, which our models ignore, might start to form in planets with an equilibrium temperature similar to that of HAT-P-67b (e.g., Lothringer et al. 2022). These clouds would lead to muted spectral features. Atmospheric temperature and metallicity can also influence the depths of the spectral lines, as they affect the extent of the atmosphere and the optical depth, respectively.

Nonetheless, the fact that we only underestimate the line depths of the ionic species (Ca ii) suggests that the lower-density atmosphere of HAT-P-67b may be more ionized than the models predict. In Figure 4, we study the effect of a higher degree of ionization on the transmission spectrum of HAT-P-67b.

Zoom In Zoom Out Reset image size

Figure 4 shows the transmission spectrum near the Ca ii infrared triplet (IRT), obtained following the steps in Section 3.5. Overlaid are model templates calculated at various temperatures with the goal of simulating the effects of different ionization levels. The model template used in our analysis, calculated at the equilibrium temperature of the planet (1903 K), underestimates the Ca ii line depths by several-fold, in agreement with the α value we found for this species. However, hotter (and therefore more ionized) models can reproduce the observed Ca ii IRT line depths, demonstrating that the upper atmosphere of HAT-P-67b may be more ionized than we predict.

Dividing the values of α by their uncertainties provides us with an estimate of the detection significance of each species (e.g., Gibson et al. 2020; Nugroho et al. 2020). The results are shown in Table 1. We find significances of 13.2 σ (Ca ii) and 4.6 σ (Na i). Different metrics are used in high-resolution spectroscopy studies to estimate the significance of a detection, including S/N from cross-correlation maps (e.g., Brogi et al. 2012), Welch's T-tests that compare the distribution of cross-correlation values in and out of the planet's trail (e.g., Brogi et al. 2012; Birkby et al. 2013), and, more recently, likelihood maps (e.g., Brogi & Line 2019; Gibson et al. 2020). Works like Cabot et al. (2019) and Spring et al. (2022) have discussed the robustness of some of these metrics. Here, we find that the significances derived from our likelihood analysis are consistent with the S/N estimates from the cross-correlation maps.

We also found a strong peak near the expected velocity of the planet in the K_p –Δv_sys maps of two additional species, Fe i and Ni i, as shown in Figure 5.

Zoom In Zoom Out Reset image size

However, the location of these peaks indicates that they are likely artifacts of the RM and CLV correction rather than true atmospheric signals. HAT-P-67b has a radial velocity semiamplitude of K_p ∼ 147 km s⁻¹, but the Fe i and Ni i signals peak at K_p ≈ 202 ± 5 km s⁻¹ (see Table 2), which coincides with the velocities at which we inject the modeled RM and CLV spectra in Section 3.1.2. Additionally, we re-ran our analysis without implementing the RM and CLV correction, and we found that the Fe i and Ni i peaks vanished. Therefore, we conclude that the Fe i and Ni i peaks are spurious signals caused by an overcorrection of the RM and CLV effects.

Visual inspection of the transmission spectra in the last row of Figure 2 reveals strong absorption near Hα, as well as strong emission and absorption near the He triplet. While the Hα signal is consistent with what we might expect from an extended hydrogen envelope around HAT-P-67b, the emission feature in the transmission spectra near the He triplet could be a manifestation of stellar variability (e.g., Palle et al. 2020b; Feinstein et al. 2021). Figure 6 shows the average, in the planet's rest frame, of all in-transit transmission spectra in the regions near Hα and the He triplet. We measure an average flux deficit at the ∼3.8% and ∼4.5% level, respectively, redshifted from the predicted line position in both cases, while the blue wing of the He feature shows an emission bump at the ∼2% level. If the absorption is planetary, these line depths would correspond to an equivalent height of the atmosphere for each species of δ R_{p,H I} ∼ 1.5R_p and δ R_{p,He I} ∼ 1.7R_p . Normalizing these values by the scale height of HAT-P-67b, we obtain δ R_{p,H I}/H ∼ 65 and δ R_{p,He I}/H ∼ 73, comparable to what Czesla et al. (2022) found for HAT-P-32b.

Zoom In Zoom Out Reset image size

We calculated transmission light curves to study the time dependence of the hydrogen and helium signals. We took an average of the relative flux of the transmission spectra within narrow passbands centered on Hα (6564.61 Å) and the He triplet (10833.22 Å) in the planet rest frame. During the averaging, we weighted each pixel by its width. The bandpasses of Hα and the He triplet had widths of 2 and 3.3 Å, respectively, such that they corresponded to equal widths in velocity space. Figure 7 presents the light curves of the Hα and He triplet signals. These light curves show that there is strong absorption before the start of the optical transit, and they justify our choice of only using the post-transit spectra as a baseline to remove the stellar features (see Section 3.5).

Zoom In Zoom Out Reset image size

We offer three possible scenarios to explain the He signal seen in the transmission spectrum of HAT-P-67b:

Scenario 1: The signal is of planetary origin, and all "out of transit" spectra are contaminated by planetary absorption. As described in Section 2, only four out of the 37 spectra in our data set were taken outside the optical transit: two before the point of first contact, and two after the point of fourth contact. In total, we have ∼20 minutes of pre-transit observations, and ∼20 minutes of post-transit data. The transmission light curves show strong absorption in the pre-transit spectra. However, the presence of planetary absorption also in the post-transit spectra (i.e., the baseline we used to remove the stellar signal) could result in a spurious emission signal. A helium outflow with a radius of R_env ∼ 2.2R_p ∼ 4.5R_J would lead to a transit in the He 10833 Å triplet that starts 20 minutes earlier and finishes 20 minutes later.

Planetary helium absorption outside the optical transit due to an outflowing envelope has been reported in other planets. For example, Nortmann et al. (2018) found post-transit absorption in the Saturn-mass planet WASP-69b that extends ∼22 minutes beyond the end of the optical transit, which they interpret as helium escape from a comet-like envelope that trails behind the planet. Similarly, Spake et al. (2021) and Zhang et al. (2023) detected significant post-egress absorption in the warm gas giant WASP-107b and the mini-Neptune TOI-2076b, respectively. Meanwhile, Czesla et al. (2022) reported redshifted absorption before the optical transit of HAT-P-32b—a gas giant that, like HAT-P-67b, orbits an F-type star on a close-in orbit. Zhang et al. (2023) also reported pre-transit absorption in the young mini-Neptune TOI 1430.01, although they do not rule out stellar variability as the source of the signal. An up-orbit stream (e.g., Lai et al. 2010; Carroll-Nellenback et al. 2017; McCann et al. 2019) can plausibly explain pre-transit absorption. In an up-orbit stream, material flows toward the star through a nozzle at the L2 point, which leads to (redshifted) pre-ingress absorption that is stronger than (blueshifted) post-egress absorption, a behavior consistent with what we see in our data. MacLeod & Oklopčić (2022) found that absorption both before and after the optical transit may be explained by moderate levels of impingement of the stellar wind on the planetary outflow. They also pointed out that, in this scenario of moderate stellar wind interaction, there is a component of absorption whose velocity is constant in the stellar rest frame. Therefore, the He signal not tracing the planetary motion does not necessarily rule out the planet as the source of the absorption.

Scenario 2: The signal is due to stellar variability. Without a longer out-of-transit baseline, we cannot rule out stellar variability as the source of the changes seen in the He line. In Figure 8 we plotted the evolution of the stellar spectrum in the region of the Si line (near 10830 Å) and the He triplet during the observations, with time going upwards. This is similar to the second panel of Figure 2, but it allows us to better visualize any changes in the stellar line profiles. Additionally, in Figure 8 we median-combined the host star spectra in trios (in pairs in the case of the out-of-transit spectra) to decrease the noise in the plot. Figure 8 shows how the median of the two pre-transit spectra (lowest line) differs greatly from the median of the two post-transit spectra (highest line). Before the transit, the stellar Si line and the He triplet form a "W" shape, while after the transit, the He triplet is significantly deformed. This means that if the signal is of stellar origin, it cannot be due to the contrast effect (e.g., Cauley et al. 2018), whereby a planet transits a heterogeneous stellar disk, with starspots and faculae, leading to spurious features in the transmission spectrum. Instead, the signal must be caused by changes in the disk-integrated stellar spectrum. The star is a fast rotator, so these changes could in principle be induced by the rotation of the stellar surface over the course of the observations.

Zoom In Zoom Out Reset image size

The change of the stellar He triplet profile outside the optical transit also rules out that the He signal we observe in the transmission spectra is due to uncorrected RM and/or CLV effects. Additionally, the contribution from these effects is well below the level of variability we measure in the data. We note that the methodology we applied to correct for these effects is not applicable to the He triplet, as photospheric models do not include chromospheric lines, and the limb-angle dependence of this line is unknown (e.g., Czesla et al. 2022). However, the contamination at the He triplet should be lower than that at Hα due to the shallower stellar line.

Figure 8 also shows how it is not only the core of the He line that changes during the observations but also the wings. In the stellar rest frame, the variability of the stellar spectra spans ∼3.25 Å, which at these wavelengths corresponds to ∼90 km s⁻¹. This value is larger than the change in radial velocity of the planet during the observations (∼60 km s⁻¹) and larger than the FWHM of observed planetary helium outflows (∼30 km s⁻¹; e.g., Zhang et al. 2023).

The Ca ii infrared triplet is a good diagnostic of chromospheric activity (e.g., Linsky et al. 1979). In our data, the Ca ii IRT lines are dominated by planetary absorption rather than showing indications of stellar activity. We know that the absorption is planetary because the Ca ii cross-correlation map, mainly driven by the deep Ca ii IRT lines in the model template, has a strong peak at precisely the radial velocity semiamplitude of the planet (Figure 3). In addition, we do not see any other significant features in the Ca ii cross-correlation map that might be caused by stellar variability.

Scenario 3: The signal is due to a combination of planetary absorption and stellar variability. The two previous scenarios are not mutually exclusive, and the He signal we see in the data could result from the combination of both.

Given that we are not able to confidently rule out stellar variability as the source of the He signal, we must also question the nature of the Hα signal. Guilluy et al. (2020) found an anticorrelation between the night-to-night variations of the out-of-transit fluxes of HD 189733 at the core of the Hα and He lines. While we find the opposite effect in HAT-P-67 (i.e., the depths of both stellar lines decrease during the observations), we do not rule out stellar variability as the source of the Hα signal, as this line is also sensitive to changes in the chromosphere. However, we can confidently claim that the source of both Hα and He signals must be astrophysical. The width and amplitude of the variability and the lack of changes in neighboring lines rule out data reduction artifacts such as the correction of telluric lines, cosmic rays, and bad pixels, or the continuum normalization.

If we assume that the He signal is caused by a planetary outflow, we can fit the transmission spectrum in Figure 6 with a one-dimensional, isothermal, H+He Parker wind model to constrain the mass-loss rate, the outflow temperature, and its line-of-sight bulk velocity. To this end, we use p-winds (Dos Santos et al. 2022), a Python implementation of the model described by Oklopčić & Hirata (2018) and Lampón et al. (2020). For our model, we assume no limb darkening and a 90% H and 10% He atmosphere. There are no available measurements of the X-ray or EUV flux of HAT-P-67, so we use the monochromatic fluxes relative to the He and H ionization edges of HAT-P-32 (Czesla et al. 2022), a star whose effective temperature (T_eff) closely resembles that of HAT-P-67. For the monochromatic flux relative to the He triplet ionization edge (911–2593 Å), we integrate the flux of a photospheric model (Castelli & Kurucz 2003) with the same T_eff, metallicity (Fe/H) and surface gravity ($\mathrm{log}g$) as HAT-P-67. In the fits, we used the spectra in a 10 Å region centered on the He feature, without masking the bump near 10832 Å. We retrieve a mass-loss rate of $\dot{M}={1.083}_{-0.058}^{+0.061}\times {10}^{13}\,{\rm{g}}\,{{\rm{s}}}^{-1}$, an outflow temperature of $T={9890}_{-340}^{+350}\,{\rm{K}}$, and an outflow line-of-sight bulk velocity of v = 13.32 ± 0.43 km s⁻¹. The posterior distributions are shown in Figure 9, and the retrieved model transmission spectrum in Figure 6.

Zoom In Zoom Out Reset image size

Assuming that the only source of heat powering the outflow is photoionization of hydrogen and helium, we can estimate the maximum mass-loss rate for our retrieved isothermal Parker wind model. Based on the calculations presented in Vissapragada et al. (2022) and implemented in p-winds, we derive a maximum mass-loss rate an order of magnitude larger than the retrieved mass-loss rate. This means that our Parker wind model is energetically self-consistent, and photoionization alone can provide enough energy to power it.

In Figure 6 we also include a Gaussian fit to the Hα signal to compare its velocity shift to that of the modeled helium outflow. From the Gaussian toy model, we obtain a Doppler shift of v = 12.36 ± 0.56 km s⁻¹, similar to the offset measured in the He signal (v = 13.32 ± 0.43 km s⁻¹). We must remember though that if the signals are of planetary origin, we are likely using a contaminated out-of-transit baseline. This means that future observations with uncontaminated out-of-transit baselines may find a different line profile, meaning that the retrieved mass-loss rate, outflow temperature, and outflow velocity will change.